Phosphorene, a monolayer of phosphorus arranged in hinge-like hexagonal structures, has attracted increasing attention owing to its intriguing thermal,^[1] electronical,^[2–4] and optical properties.^[5] The successful exfoliation of phosphorene from bulk black phosphorus has spurred tremendous research interest into this two-dimensional (2D) material since 2014.^[6–8] Its unique structure and intriguing anisotropic properties are rarely found in other 2D materials. Being chemically and electronically active, its potential applications such as photodetectors,^[9–12] field-effect transistors,^[13–22] rechargeable batteries,^[23–28] and thermoelectric devices^[29–31] have also drawn growing research interest. Due to its ultrahigh surface/volume ratio and high chemical activity, pristine phosphorene has been shown to be a promising candidate for gas sensing.^[32] It has been revealed that by adding anionic dopants such as O, C, and S atoms, a large cohesive energy and highly dispersive electronic states can be formed in doped phosphorene.^[33] The band gap decreases by pushing down the conduction band, suggesting that the optical and electronic properties of the doped phosphorene can be tailored. Particularly, the metal doped phosphorene exhibits significantly enhanced chemical activity compared with pristine phosphorene, suggesting its great potentials in chemical applications such as materials growth, catalysis, gas sensing and storage. Moreover, due to the peculiar puckered structure of phosphorene, a negative Poissonʼs ratio is observed in the out-of-plane direction under uniaxial deformation in the direction parallel to the pucker. This phenomenon originates from the in-plane lattice puckered structure, together with coupled hinge-like bonding configurations, and has been confirmed by both numerical^[34] and experimental^[35] studies.

Besides pristine phosphorene, novel phosphorene-based heterostructures have also shown great potentials for electronic and optoelectronic applications. It has been found that many desirable electronic characteristics of phosphorene such as direct band gap and linear dichroism can be preserved in heterostructures such as phosphorene/graphene and phosphorene/hexagonal boron nitride (h-BN) bilayers.^[36] Graphene and h-BN sheets can be used as effective capping layers to protect phosphorene from structural and chemical degradations. Meanwhile, both sheets can also act as active layers to tune the carrier dynamics and optical properties of phosphorene. A recent study based on first-principles calculations by Gao et al.^[37] has shown that the stability of phosphorene nanoflake is strongly dependent on the coupling strength between phosphorene and its substrate. For example, a strong interaction of 0.75 eV/P atom with copper substrate can break up the phosphorene structures, while a weak interaction of 0.063 eV/P atom with h-BN cannot stabilize phosphorene. However, it is reported that a moderate interaction strength of ∼0.35 eV/P atom is able to keep the 2D characteristics of phosphorene nanoflake on a realistic time scale. Therefore, in the study of phosphorene-based heterostructures, it is important to seek the optimum material and coupling strength to avoid the edge wrapping and reconstruction of phosphorene structures.

In this article, the state-of-the-art progress on the numerical studies of thermal transport properties of 2D phosphorene and phosphorene-based materials are reviewed, aiming to summarize the up-to-date results for this emerging material in the field of energy transport research. This review is organized as follows. First, the fundamental mechanisms and basic principles of different numerical methods for characterizing thermal conductivity (κ) and interfacial thermal resistance (R) are introduced. Next, detailed descriptions of the in-plane and out-of-plane thermal transport in phosphorene and phosphorene-based heterostructure materials are classified and highlighted, followed by a discussion of their potential applications in electronic devices such as photodetectors, field-effect transistors, lithium ion batteries, sodium ion batteries, and thermoelectric devices. Lastly, a summary of the numerical research in the field of thermal transport in 2D phosphorene is highlighted along with our perspective on its great potentials and opportunities in this research area.

Boltzmann transport equation (BTE) associated with the first-principles calculations have long been used to calculate the thermal conductivities of crystalline materials. The first-principles lattice thermal conductivity calculations for monolayer phosphorene start with the computation of interatomic force constants (IFCs) based on the density functional theory (DFT). Next, the dynamical matrices are used to solve the BTE with relaxation time approximation (RTA),^[38] in which the lattice thermal conductivity is expressed as

Another numerical approach for thermal transport property characterization is the non-equilibrium Green function (NEGF) method, which can exactly deal with ballistic thermal transport and gives the maximum thermal conductance of a material.^[49–51] When the system size of a nanomaterial is smaller than its phonon mean free path, the thermal transport can be regarded as ballistic, which can be described by the Landauer formula

Molecular dynamics (MD) simulation is another powerful tool to treat thermal transport problems at micro/nanoscale. It intrinsically includes full anharmonicity in atomic interactions and does not make any assumptions on the thermodynamic limit. Since the description of atomic trajectories is achieved by numerically solving Newtonʼs equations of motion, the MD method can deal with thermal transport problems in systems containing millions of atoms. Therefore, MD simulations have been used in both thermal conductivity and interfacial thermal resistance research.^[52–59] Tremendous efforts have been devoted to developing empirical interatomic potential (EIP) fields that can be adopted in MD simulation. One common strategy to develop an EIP is to first obtain the materials properties, e.g., crystal structure, cohesive energy, or phonon dispersion, from either first-principles calculations or experimental measurements, and then parameterize the potential by best fitting those properties.^[60,61] Specific to the MD simulations in phosphorene, several types of EIPs have been developed to describe the P–P interactions. Jiang et al.^[62] parameterized the Stillinger–Weber (SW) potential to describe the single layer phosphorene, which is expressed as^[63]

The methods to characterize the thermal conductivity of phosphorene using MD simulations can be categorized into three groups, i.e., the steady-state equilibrium method, steady-state non-equilibrium method, and transient method. The equilibrium molecular dynamics (EMD) method is also referred to as the Green–Kubo method, which is given by the expression of heat current autocorrelation function (HCAF)^[69]

Aside from the in-plane thermal conductivity, the cross-plane thermal contact resistance is also an important attribute often considered in thermal interface materials. As aforementioned, the NEMD method has been traditionally used to characterize the thermal conductivity of micro/nanostructures. Besides, the NEMD method has also been widely applied to calculate the interfacial thermal resistance between adjacent materials. After the hybrid system reaches thermal equilibrium, a heat flux is imposed on the system which will flow across the contact boundary. The temperature jump at the interface can be used to calculate the thermal contact resistance based on Fourierʼs law of heat conduction. However, if the 2D phosphorene is directly placed on the substrate, the NEMD method cannot be applied to calculate the interfacial thermal resistance between them. To build a steady-state temperature gradient in the substrate and across the interface, the heat reservoir needs to be placed directly on this phosphorene layer. In MD simulations, the temperatures are calculated according to the energy equipartition theorem

To remedy this problem, a transient method can be used to calculate the interfacial thermal resistance between supported phosphorene and substrate without measuring and controlling the temperature of phosphorene simultaneously. Basic principles of this method are illustrated as follows. First, temperature controls are applied to different components of the heterostructure to reach thermal equilibrium at steady-state. Then, an ultrafast thermal impulse is imposed on the phosphorene to elevate its temperature to a much higher value. In the meantime, the surface temperature of the substrate material can be considered as unchanged since the excitation time is very short. Next, in the thermal relaxation process, thermal energies in phosphorene will be gradually transferred to the substrate through heat conduction until the two systems have reached thermal equilibrium again. Since the only thermal pathway for heat dissipation is through heat conduction from phosphorene to substrate, the total energy variations of phosphorene can be correlated with the temperature evolutions in phosphorene and substrate surface, which is expressed as

The puckered atomic structures of monolayer phosphorene are shown in Figs. 1(a) and 1(b). Since thermal conductivity of most semiconductors is mediated mainly by lattice vibrations, the calculated κ is directly related to the crystals’ phonon dispersion relations.^[85–87] Similar to other 2D materials such as graphene and silicene,^[88–90] monolayer phosphorene also has a quadratic flexural acoustic (ZA) phonon branch near the Γ point, as shown in Fig. 1(c). The ZA phonons in phosphorene vibrate perpendicularly to the basal plane, which is similar to that in graphene. On the other hand, the ZA phonons in silicene do not vibrate solely in the out-of-plane direction.^[89] It can be observed from Fig. 1(c) that the longitudinal acoustic (LA) phonon branch is asymmetric around the Γ point. Based on the slopes of the LA phonon branches near the Γ point, the group velocities along the zigzag Γ–X and the armchair Γ–Y chiralities are predicted to be 7.83 km/s and 4.01 km/s, respectively.^[91,92] In a different study, Zhu et al.^[39] calculated the group velocities of the LA phonon modes along zigzag and armchair directions to be 8.6 km/s and 4.5 km/s. The experimentally measured group velocities of bulk phosphorus along the zigzag and armchair directions equal 9.6 km/s and 4.6 km/s, respectively.^[93]

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. (color online) (a) Top view of a 3 × 3 phosphorene supercell, the zigzag and armchair chiralities are along the x and y directions, respectively. (b) Enlarged single unit cell of phosphorene. The top and bottom phosphorus atoms are marked separately as P1 and P2. Reproduced from Ref. [65] with permission from AIP Publishing. (c) Phonon dispersion relation of monolayer phosphorene along the path passing through high-symmetry k-points. The acoustic and optical phonon branches are denoted in different colors. The upper inset denotes the high-symmetry k-points in the Brillouin zone and the lower inset depicts the vibration direction of flexural phonons. Reproduced from Ref. [91] with permission from the Royal Society of Chemistry.

The anisotropic thermal conductivity of phosphorene can be attributed to the discrepancies in the group velocities, which is possibly caused by the different atomic arrangements in the hinge-like structures of phosphorene.^[94] The anisotropic thermal conductivity has also been reported in other hinge-like structures such as SnSe,^[95] SnS, and arsenene.^[96] Aside from the pronounced spatial-anisotropy of thermal conductivity in phosphorene, first-principles calculations also reveal that phosphorene possesses an anisotropic electrical conductance, i.e., the electronic conductivity along the armchair direction is much higher than that along the zigzag direction.^[97] The prominent electrical and thermal conducting directions are orthogonal to each other, enhancing the thermoelectric figure of merit. Particularly, the prominent electronic conductivity combined with lowered thermal conductivity in the armchair direction indicates that the phosphorene nanoribbon might be a promising thermoelectric candidate.

The anisotropic thermal transport in phosphorene has been quantitatively assessed by several first-principles studies. By solving the BTE based on first-principles calculations, Qin et al.^[91] calculated the anisotropic thermal conductivity of monolayer phosphorene at different temperatures. The predicted κ at 300 K in zigzag and armchair directions are and , respectively, showing an obvious anisotropy along different directions. It is reported that in all temperature ranges, the thermal conductivity in the zigzag direction is generally twice the value of that in the armchair direction. Besides, the calculated κ has a perfect inverse relationship with temperature when T is higher than the Debye temperature (278.66 K), indicating that Umklapp scattering is dominant in this temperature range. The calculated lattice thermal conductivities of armchair and zigzag phosphorene at different temperatures are shown in Fig. 2. Also using first-principles calculations combined with the non-equilibrium Green function (NEGF) method, Ong et al.^[98] investigated the ballistic thermal transport in monolayer phosphorene. A similar conclusion has been drawn that the thermal conductance of phosphorene along the zigzag direction is about 40% higher than that along the armchair direction. Besides, owing to the peculiar atomic structures in phosphorene, they reported an anomalous effect of uniaxial strain on the predicted thermal conductance. In stark contrast, it is discovered that the thermal conductance in the zigzag direction can be enhanced when a parallel in-plane strain is applied, but decreases when a perpendicular in-plane strain is applied. On the other hand, the thermal conductance in the armchair direction always decreases with either parallel or perpendicular strains. The changes from the contributions of acoustic and optical phonons account for the peculiar dependence of the thermal conductance with uniaxial strains. The strain engineering of the thermal conductivity of phosphorene has also been confirmed by Cai et al.^[99] using detailed first-principles calculations. By examining the whole 2D Brillouin zone, an extraordinary phononic anisotropy in phosphorene was demonstrated, which depends sensitively on the magnitude and type of the applied strain. Using DFT calculations combined with BTE, Zhu et al.^[39] revealed a compelling coexistence of size-dependent and size-independent thermal conductivities in phosphorene. The calculated thermal conductivities in the zigzag and armchair chirality equal and at temperature 300 K. They discovered that intriguingly, the thermal conductivity along the armchair direction of phosphorene is almost a constant ( ) with respect to the q-mesh size, while the thermal conductivity along the zigzag direction steadily increases with the number of q points. The close to zero group velocities for TA phonons in the armchair direction and the periodically corrugated structure stacking along the armchair direction are used to explain their calculation results. Also using first-principles calculations combined with BTE, Liu et al.^[100] systematically calculated the anisotropic thermal conductivity in phosphorene in a complete set of crystal chirality, i.e., transport direction, ranging from 0.0° (zigzag) to 90.0° (armchair). It was reported that the intrinsic thermal conductivity is a smooth monotonic decreasing function of the crystal chirality, which exhibits a sinusoidal behavior bounded by the two terminated values (0.0°) and (90.0°). Besides, they illustrated that the optical modes in phosphorene have unusually large contributions to heat transfer, which account for almost 30% of the total thermal conductivity. The large contributions from optical phonons are attributed to the comparable group velocities and relaxation time with the acoustic phonons. The orientation-dependent thermal conductivity of phosphorene is shown in Fig. 3.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Lattice thermal conductivity of phosphorene along zigzag and armchair directions with temperature ranging from 200 K to 800 K. The κ fitted by the inverse relationship with temperature ( ) is plotted with dot lines. Reproduced from Ref. [91] with permission from the Royal Society of Chemistry.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Orientation-dependent thermal conductivity of phosphorene ranging from 0.0° (zigzag) to 90.0° (armchair). The gray dot line is fitted by a cosine function with a period of 2π. Reproduced from Ref. [100] with permission from the Royal Society of Chemistry.

Aside from the aforementioned studies using first-principles calculations, MD simulations have also been performed to calculate the thermal conductivity of phosphorene. A recent study by Hong et al.^[101] also confirmed that the thermal conductivity of phosphorene in the zigzag direction is much higher than that in the armchair direction. By using the NEMD method, the extrapolated thermal conductivities for infinite length phosphorene in zigzag and armchair directions equal and , respectively. A typical setup of NEMD simulation and the fitting method on temperature profiles are shown in Figs. 4(a) and 4(b). They attributed the strong anisotropic thermal transport to the distinct atomic structures at altered chiral directions and direction-dependent group velocities. Another MD simulation study by Xu et al.^[65] using the EMD method predicted a much larger anisotropy ratio of five. The calculated thermal conductivities along zigzag and armchair directions of phosphorene are and , respectively. A weak anisotropy of phonon lifetimes is revealed by lattice dynamics calculations. Besides, the effective phonon mean free paths at zigzag and armchair directions are extracted using the relaxation time approximation (RTA), which equal 141.4 nm and 43.4 nm, respectively. In-plane thermal conductivity of phosphorene phononic structures has also been investigated by using the EMD method.^[102] It is reported that when the porosity in the phononic crystal reaches 36%, the thermal conductivity in both zigzag and armchair directions is reduced to below 5% of their original values for pristine structures. The remarkable reduction of thermal conductivity could lead to promising applications of phosphorene phononic structures in thermoelectric devices. The strong reductions of thermal conductivities are originated from the depression of phonon group velocities induced by Brillouin zone folding and the diminished phonon lifetimes.

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. (color online) (a) Illustration of the NEMD process. Black atoms at the boundaries of the system are fixed in the position. Red and blue areas are denoted as heat bath and heat sink, respectively. (b) Temperature distribution along the heat flux direction in the armchair direction. The red solid line denotes the linear fitting results. Reproduced from Ref. [101] with permission from the Royal Society of Chemistry.

Although the existing numerical studies have mostly agreed that the thermal conductivity in phosphorene is highly anisotropic, however, some researchers have posed questions on this open topic. Using NEGF and first-principles method, Jiang et al.^[103] investigated the thermal conductance for monolayer phosphorene in the ballistic transport regime. It is found that the anisotropy in the thermal conduction is very weak, with only 4% differences between armchair and zigzag directions. Detailed phonon dispersion calculations disclose that the ZA mode has lower group velocity in the direction perpendicular to the pucker, while the LA mode has higher group velocity in the perpendicular direction. This unusual phenomenon is attributed to the highly anisotropic Poisson ratio in phosphorene. As a result, the competition between these two opposite effects leads to the weak anisotropic thermal conductance for phosphorene.

Even though there are discrepancies in the calculated thermal conductivity of phosphorene at different chiral directions, it is affirmative that the in-plane thermal conductivity of phosphorene is about an order of magnitude smaller than that of graphene.^[101] The significant thermal conductivity differences between phosphorene and graphene can be attributed to their different phonons behaviors in the out-of-plane direction. By examining the contributions of different phonon branches to the thermal conductivity, it is reported that the ZA phonons in phosphorene only contribute to ∼5% of the overall thermal conductivity, which is much lower than that of graphene.^[91] Recent studies have shown that for a suspended single layer graphene, the ZA phonons can contribute up to 77% at 300 K and 86% at 100 K of the thermal conductivity due to its high specific heat and long phonon scattering time.^[104] Besides, by numerically solving the phonon BTE, Lindsay et al.^[88] came up with a symmetry-based selection rule which significantly restricts anharmonic phonon–phonon scattering of the ZA phonons, and they proved that κ of graphene is dominated by the ZA phonon modes. Hong et al.^[101] compared the thermal conductivity differences between phosphorene and graphene. The overall and decomposed phonon density of states (PDOS) are calculated in both structures to elucidate their different phonon behaviors. The calculated results are shown in Fig. 5. The plotted PDOS of graphene coincides well with that calculated from first-principles,^[105–107] which illustrates that the ZA branch dominates the low frequency acoustic phonons while the in-plane LA and TA branches occupy the high frequency phonons. On the other hand, the PDOS in phosphorene is isotropic in all directions, i.e., the ZA phonons occupy the same frequency range as those of LA/TA. Moreover, it can be clearly observed that vibrational modes excited in phosphorene are severely limited compared with graphene. The active phonon frequencies of the phosphorene sheet range from 0 to 15 THz, while the lateral phonons in graphene dominate the high frequencies up to 52 THz and the ZA phonon occupies the low frequency acoustical branches. Based on the above discussions, it can be concluded that the lower thermal conductivity of phosphorene compared to that of graphene could be caused by the minor contributions from the ZA phonons modes.

Fig. 5.

Figure Option ViewDownloadNew Window

Fig. 5. (color online) (a) Overall and decomposed phonon density of states in (b) x direction, (c) y direction, and (d) z direction for phosphorene and grapheme respectively. Two phenomena are observed. First, the PDOS of graphene is anisotropic, while that of phosphorene is isotropic. Second, the active phonon frequency in phosphorene is severely limited compared to that of graphene. Reproduced from Ref. [101] with permission from the Royal Society of Chemistry.

It has been widely acknowledged that 2D phosphorene monolayers are chemically active and unstable when exposed to ambient environments.^[108] Recent studies have shown that being placed on a substrate or sandwiched in matrix materials could greatly enhance its structural and thermal stability.^[109] The thermal performance of these phosphorene-based heterostructures needs to be further investigated. The interfacial thermal transport between monolayer phosphorene and crystalline silicon substrate has been characterized by Zhang et al.^[83] using the transient pump-probe method as aforementioned. Illustration of the pump-probe technique and schematic of the hybrid system setup are shown in Figs. 6(a) and 6(b). It was revealed by MD simulations that at room temperature, the interfacial thermal resistance between phosphorene and silicon is very low and independent of the silicon thickness when the substrate thickness is larger than 3.12 nm. Effects of different external and internal modulators such as system temperature and coupling strength have also been investigated and it was reported that R decreases remarkably with increasing temperature and contact strength. The superiority of phosphorene on interfacial heat transfer was further proved by comparing its R values with its 2D counterparts graphene and silicene. Two distinct advantages of the phosphorene/silicon over graphene/silicon and silicene/silicon are discovered. Firstly, within the studied ranges of system temperature and coupling strength, the interfacial thermal resistance between phosphorene/silicon is about a quarter of that between graphene/silicon, which proves that 2D phosphorene has better heat dissipation performance in interfacial heat transfer compared with graphene. Secondly, with the increases in temperature and coupling strength, the interfacial thermal resistance between phosphorene/silicon decreases more sharply than that between silicene/silicon, indicating that phosphorene is more sensitive to ambient condition changes.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) (a) Schematic of the transient pump-probe process. (b) Atomic configuration of a phosphorene/silicon hybrid system. Reproduced from Ref. [83] with permission from ScienceDirect.

Lately, it was shown that stacking a graphene/phosphorene van der Waals (vdW) bilayer can preserve their properties in the ultimate heterostructure.^[110] The relative position of phosphoreneʼs band structure with respect to grapheneʼs can be tuned via a vertical external electric field. Moreover, by exploring the electric field dependent band structures and optical properties of the graphene/phosphorene bilayer system, Hashmi et al.^[111] demonstrated that the bilayer heterostructure can be applied to a high-speed device although the optical anisotropy in the bilayer structure for in-plane electric field polarization has disappeared. Due to the presence of the lone-pair state, monolayer phosphorene can be corrugated when in contact with common metal electrodes, which may degrade its performance. Conversely, graphene has excellent structural integrity with both metal electrodes and phosphorene due to its atomically smooth surface. Thus, graphene can serve as a perfect interfacial material between the phosphorene and metal electrodes.^[36] Owing to the above reasons, thermal transport in graphene and phosphorene composed heterostructures have attracted tremendous research interests. Extensive MD simulations have been performed by Zhang et al.^[112] to investigate the interfacial thermal transport in multilayer graphene/phosphorene/graphene heterostructures along the cross-plane direction. The hybrid system setup is shown in Fig. 7(a). Based on the simulation results, the interfacial thermal conductance at the graphene/phosphorene interface can be effectively controlled by several modulators such as vacancy defects, hydrogenation, and cross-plane strain. At the largest defect concentration of 5% under consideration, the hydrogenation and vacancy defects created in the interfacial graphene layer can bring 3.2-fold and 1.7-fold increases in interfacial thermal conductance, respectively. Effects of defect and hydrogenation on the calculated interfacial thermal conductance are shown in Fig. 7(b). This defect-induced enhancement of the interfacial thermal conductance is mainly attributed to two reasons. Firstly, the elastic modulus mismatch decreases with increasing defect level. Secondly, the interlayer interaction strength between graphene and phosphorene is enhanced. External cross-plane strains can also be used to manipulate the interfacial thermal conductance between graphene/phosphorene/graphene heterostructures. A cross-plane compressive strain of 0.05 can increase the interfacial thermal conductance by 90%, while a cross-plane tensile strain of 0.03 can reduce the interfacial thermal conductance by 30%. In another work, Hong et al.^[76] calculated the out-of-plane thermal transport properties between bilayer phosphorene and graphene using the pump-probe method. The hybrid system setup is shown in Fig. 8(a). Several modulators such as system temperature, contact pressure, vacancy defect, and hydrogenation are explored, with which significant thermal resistance reductions are observed. The maximum R reduction is predicted as 84% when the hydrogenation is applied on the near-phosphorene-side graphene surface. Other factors such as temperature, coupling strength, and fraction of single-vacancy defects have relatively weaker influences on R, which decrease R values by 57%, 70%, and 35%, respectively. Effects of temperature and interatomic coupling strength on the interfacial thermal resistance are shown in Fig. 8(b). The PDOS mismatch in graphene and phosphorene appears to be the key factor for the thermal resistance. Reductions of R at the interface can be attributed to several factors, including increased phonon population, enhanced anharmonic phonon scattering at higher temperatures, as well as strengthened coupling between lateral and out-of-plane phonons of graphene with increasing fraction of the defect or functionalization. Note also that the interfacial thermal resistance between phosphorene/graphene is smaller than that of other 2D bilayer heterostructures, which is a distinct advantage for practical applications such as high-performance thermal radiators.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) (a) Simulation model of graphene/phosphorene/graphene heterostructure in which a 16-layer phosphorene is sandwiched in an 8-layer graphene at both ends. (b) Interfacial thermal conductance as a function of defect concentration and hydrogenation. Reproduced from Ref. [112] with permission from IOP Science.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (color online) (a) Atomic configuration of the phosphorene/graphene bilayer heterostructure. (b) Dependence of the interfacial thermal resistance on the temperature and coupling strength. The predicted R decreases monotonically with temperature and the contact pressure. Reproduced from Ref. [76] with permission from the Royal Society of Chemistry.

Also using MD simulation method, Pei et al.^[113] investigated the thermal stability and thermal conductivity of phosphorene in phosphorene/graphene heterostructures. It is found that when phosphorene is encapsulated in-between two graphene sheets, its thermal stability is enhanced significantly. The thermally stable temperature of phosphorene shows an increase of 150 K. However, no increase in the thermally stable temperature is observed when phosphorene is on a graphene sheet. Moreover, it was reported that the thermal conductivity of phosphorene in phosphorene/graphene heterostructures is much higher than that of a free-standing one. A net increase of 20%–60% in thermal conductivity is obtained from their simulations. The thermal conductivity alternations for pristine phosphorene, phosphorene/graphene bilayer, and graphene/phosphorene/graphene sandwiched structures are shown in Fig. 9. Phonon spectral energy density analysis suggests that the increased thermal conductivity of phosphorene in phosphorene/graphene heterostructures is mainly due to the enhanced group velocities and extremely strong phonon coupling between phosphorene and the graphene substrate. In a similar study, Chen et al.^[114] investigated the interfacial thermal conductance in graphene/black phosphorus heterogeneous structures using the MD simulation method. Two important parameters, i.e., the critical heat power density and the maximum heat power density, are identified. It was shown that the interfacial thermal conductance can be effectively tuned in a wide range by external strains and interface defects. The compressive strain can enhance the interfacial thermal conductance by one order of magnitude, while interface defects give a two-fold increase. Recently, using the NEMD method, Zhang et al.^[115] studied the in-plane and cross-plane thermal conductivities of single- and multi-layer phosphorene. The variation in the thermal conductivity with respect to the sample length, orientation, layer number, and mechanical strain was systematically examined. From the simulation results, it was found that the phosphorene possesses a strong anisotropy and also layer-number independence of the in-plane thermal conductivity due to its puckered structure, which is distinctly different from other 2D materials. Similar to other 2D materials, the in-plane thermal conductivity is dependent on the sample length and mechanical strain. Compared with the in-plane thermal conductivity under uniaxial in-plane strain, the cross-plane thermal conductivity is more sensitive to the cross-plane strain and is thus more tunable by the mechanical strain.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (color online) Thermal conductivities in the (a) armchair and (b) zigzag directions of phosphorene for the three models, i.e., free-standing phosphorene, phosphorene on a graphene sheet, and phosphorene sandwiched between two graphene sheets. Reproduced from Ref. [113] with permission from the PCCP Owner Societies.

This paper offers an overview of recent numerical studies on the anisotropic thermal conductivity of phosphorene and the in-plane and cross-plane thermal transport in phosphorene-based materials. Due to their unique thermal and electrical characteristics, phosphorene-based materials can offer outstanding performance for various applications such as photodetectors, field-effect transistors, lithium ion batteries, sodium ion batteries, and thermoelectric devices. In spite of its great potential to be used in electronic and optoelectronic applications, challenges remain in the implementations of phosphorene in these devices. Due to its high surface area and volume ratio compared with bulk black phosphorus, its chemical and thermal stability is modified. In contrast to the 550 °C sublimation temperature of bulk phosphorus, the decomposition temperature of phosphorene is 400 °C.^[108] Besides, a recent experimental study has shown that the volume of few layer phosphorene can be increased by 200% when exposed in ambient environments due to condensation of moisture from air.^[116] Furthermore, long term exposure to ambient conditions can result in a layer-by-layer etching process in phosphorene flakes, which could transform few layer flakes into single layer phosphorene. In the time scale of a few minutes, a shift in the threshold voltage of phosphorene occurs due to the physisorbed oxygen and nitrogen. In longer time scales of hours, the p-type doping occurs from water adsorption. The electronic properties of phosphorene are quite sensitive to external turbulence due to its strong affinity of water molecules. An itinerant behavior of mono- and di-vacancies in phosphorene is observed at room temperature.^[117] The mobility of a mono-vacancy is 16 orders faster than that in graphene, and its hopping rate in the zigzag direction is ∼3 orders higher than that in the armchair direction. The exceedingly high motilities of mono-vacancy and di-vacancies can be used to explain the relatively low structural stability and high affinity to environmental molecules in phosphorene. One possible way to suppress the oxidation of phosphorene in an ambient environment is to embed it in van der Waals heterostructures and apply a vertical electrical field. For example, Gao et al.^[109] reported that phosphorene/MoSe₂ heterostructure is able to reverse the stability of physisorption and chemisorption of O₂ molecular on phosphorene. The application of a vertical electric field of −0.6 V/Å can further increase the oxidation energy barrier to 0.91 eV. However, the effectiveness of this method on prohibiting the adsorption of other molecules has not been verified. Despite the side effects caused by ambient environments, on the other hand, the nontrivial and distinct charge transfer occurring between phosphorene and physisorbed molecules makes it a promising candidate for gas sensor applications.^[118] Due to the above mentioned reasons, further studies need to be undertaken to investigate the variations of the thermal transport in phosphorene with the existence of ambient molecular species such as H₂O, CO, H₂, NH₃, NO, NO₂, and O₂. Also, it is important to further explore phosphorene-based heterostructures that can increase its chemical and structural stability in an ambient environment while preserving or enhancing its thermal performance. Since numerical approaches such as first-principles and molecular dynamics are more cost-effective and efficient compared to experimental studies in the field of heterostructure exploration and property probing, more efforts need to be made through numerical studies to obtain the properties of phosphorene-based compositions and to promote the usage of 2D phosphorene in future electronic applications.